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ABSTRACT 


A quantitative evaluation of the limitations of the radiation boundary 
elements in the finite element code ATILA [Ref. 1] has been performed. Five 
three-dimensional models were employed, each representing a rigid spherical 
solid surrounded by water. Monopolar, dipolar and quadrupolar incident 
spherical waves were introduced and the corresponding scattered waves were 
computed using the ATILA code and an exact analytical solution. 

The dimensionless parameters that characterize the problem are ka, KL, 
and kR where k is the wavenumber of sound in water, a is the radius of the 
scatterer, R is the outer fluid mesh radius, and L is the thickness of the fluid 
layer. The range of values investigated were KR=1.5, 2.5, 4.0, ka=0.5, 1.0, 2.0 
and kL=0.5, 1.0. 

For axially symmetric incident fields, the maximum normalized errors 
occurred at the poles and were 9%, 12%, and 6%, respectively. Furthermore, 
the errors for monopolar and dipolar incident fields were strongly influenced by 
the location of the radiation boundary (kR), less so by the scatterer’s radius (ka); 
specifically, the error decreases with increasing kR and/or ka. The errors for 
quadrupolar incident fields do not exhibit any significant dependence on kR or 
Ka. The errors for all the axially symmetric incident fields were not affected by 
variations of the element’s size (kL). For non-axially symmetric incident fields, 
the maximum deviation occurred at the equatorial points and was less than 
9.5%. 

Further investigation using a two-dimensional model is proposed in order 
to determine the range of values of ka, kL, and kR which will result in negligibly 


small errors. 
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l. INTRODUCTION 


The investigation described in this thesis is part of an Ongoing research 
project in the numerical modeling of arbitrary densely packed, random volumetric 
active sonar arrays [Ref. 1]. The technique employed is an extension of the so- 
called T-matrix method, which has been applied to other scattering problems 
{Ref. 2]. This method rigorously accounts for multiple scattering to all orders. To 
apply the T-matrix method to the problem of active sonar array performance 
prediction, it is necessary to compute the radiating and scattering properties of a 
single array element in a free field environment. To accomplish this for a real 
transducer, the ATILA [Ref. 3] finite element code is employed. 

ATILA is a French finite element code especially developed for the 
analysis of underwater acoustic transducers. This code was written by engineers 
at the Institut Superieur d'Electronique du Nord (ISEN), Lille, France and is in 
use by U.S. scientists within and outside the Navy working on U.S. Navy- 
sponsored research. 

Computation of the radiating and scattering properties of a sonar 
transducer using ATILA involves building a finite-element model representing the 
transducer and surrounding it by a finite-element mesh representing water, which 
is terminated by so-called “radiation boundary elements”. The radiation 
boundary elements are intended to “absorb” all incident acoustic radiation. In 
practice, they perform this function less than perfectly. This has a direct effect 
on the computation of a transducer’s radiating and scattering properties. 

The present investigation focuses on the evaluation of ATILA's radiation 
boundary elements. For this purpose, a three dimensional spherical fluid mesh, 
surrounding a spherical rigid body and terminated by radiation boundary 
elements is studied. An incident spherical wave consisting of a single 
component of a multipolar acoustic field strikes the sphere, and the resulting 


scattered wave field is calculated using ATILA. The amounts of acoustic field in 
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the intended outgoing component and in other multipolar components are then 
analyzed. The influence on the results of the fluid mesh element size, fluid mesh 
inner radius, and fluid mesh outer radius relative to the acoustic wavelength are 
investigated and quantitative guidelines developed in order to minimize the effect 
of imperfect absorption. 

The remainder of this thesis is divided into four chapters. Chapter II 
describes the theory involved in the finite element analysis, the possible types of 
analyses that can be performed, particularly the harmonic analysis of a radiating 
spherical structure excited by an incident spherical wave. 

Chapter Ill describes the three dimensional spherical models which were 
employed. Chapter IV presents and discusses the scattering results for different 
incident multipolar components and the influence of element size and mesh 
inner and outer radius. Chapter V concludes the thesis. Appendix A contains 
the C program used to analytically calculate the pressure scattered from the 
spherical boundary for pressure release and rigid surface boundary conditions. 
Appendix B presents the FORTRAN code used in ATILA for the generation of 


the spherical incident pressure wave. 


ll. THEORY 


A. GENERAL PRESENTATION OF ATILA VERSION 5.03 

The finite element method is a technique that provides numerical 
solutions for boundary value problems and field calculations. [Refs. 4,5,6,7,8] 
ATILA is a finite element code developed specifically for the analysis of sonar 
transducers. 

The ATILA code is able to perform: 


alk elastic or piezoelectric structures modeling, 
Zz magnetostrictive structures modeling, 
3. periodic structures modeling, 


following a general formulation shown in Figure 1 (after [Ref. 3]), to model an 


underwater acoustic transducer. 


Open Exterior Fluid Domain 
Closed Interior 


Fluid Domain Fluid - Structure 
Interfaces 


IN 


Non Reflecting 
Boundary S 


Elastic 


Domain Active Domain 





Figure 1. ATILA general formulation. After Ref. [3] 


The ATILA code is based on the separation of the physical problem under 
consideration into a discrete number of elements which are in turn described by 
their nodes, in a given order. For each node there are a number of degrees of 
freedom (d.o.f.) that can be specified using certain boundary conditions. The 
elements, nodes, and the specific node-numbering order, is referred to as the 
topology of the problem. 

An ATILA job organization is carried out in several steps, as follows: 

Model definition. 
Mesh generation. 
Data file preparation. 


Running a job. 


oS SS = 


Result file postprocessing. 

First of all, the type of analysis to be performed has to be specified, for 
example, harmonic analysis of a solid structure excited by an impinging wave. 
The type of elements that are required to describe the fluid and structure domain 
are then chosen. ATILA includes a library of 46 different elements to model 
composite, piezoelectric, fluid, magnetostrictive, coupling FEM-BEM, interface, 
and radiation dampers. 

The mesh generation procedure includes the assignment of coordinates 
for each node and the node-numbering order for each element. Throughout this 
procedure the whole physical problem is discretized into elements which allow 
the representation and modeling of many different geometrical shapes and lines, 
as, for example, PRIS15F, a fifteen-node isoparametric triangular base prism 
used to model homogenous fluid media, or TRIAO6R, a six node isoparametric 
triangular element to prescribe a monopolar, dipolar, or multipolar radiation 
condition. 

One of the facilities of the ATILA code is the MOSAIQUE preprocessor. 
This routine enables the use of super-elements and generates all the necessary 


elements and node data for ATILA. 


The data file includes all the neccessary input data such as the type of 
analysis, the material properties, the node coordinates, the elements, the 
boundary conditions, the loading data, and possibly the plane wave data, to carry 
out the analysis. 

Running an ATILA job involves calling up the subroutines to compute 
elementary matrices, solve equations and display the results. 

The available results file, if desired, can be postprocessed to create 
graphic displays of the structure, contours of constant value for potentials, 
pressures, and displacements. A simple flow chart of an ATILA job is shown in 
Figure 2 (after [Ref. 3]). | 


1.MODEL DEFINITION 


2. MESH GENERATION 


3. DATA FILE PREPARATION 


4. RUNNNING A JOB 


5. RESULT FILE / POSTPROCESSING 





Figure 2. Flow chart of an ATILA job. After Ref.[3] 


B. ATILA FINITE ELEMENT COMPUTATION FOR RADIATING AND 

SCATTERING PROBLEMS 

1. General Elastic or Piezoelectric Structures Modeling 

A large number of analyses can be performed. These are based upon the 
complete set of equations of elasticity in the structure, the Helmholtz equation in 
the fluid, and Poisson's equation in the elastic or piezoelectric material. 
Appropriate radiation boundary conditions are applied on the surface which 
surrounds the fluid domain. 

The unknown quantities are the nodal values of the displacement field U 
in the whole structure, the electric potential ® in the piezoelectric material, and 


the pressure P in the fluid. The system of equations is written in matrix form: 


((K,,J-@2(Ml) {Ky}  -{L] Ul fF 
Ka in [0] ° =| -q 
—pe%o2[L]" —-(0]" {HJ -@2(M,J) | LP} Lec* 


where: 
- vector of the nodal values of the components of the displacement field. 
: vector of the nodal values of the electric potential. 


: vector of the nodal values of the pressure field. 


Im 1D Je IC 


- vector of the nodal values of the applied forces. 


- vector of the nodal values of the electric charges. 


to} 


w: vector of the nodal values of the integrated normal derivative of the 
pressure on the surface boundary S. 


(K Stiffness matrix. 


ip ; 
[K,,]; piezoelectric matrix. 
[K,,]: dielectric matrix. 


[M]: consistent mass matrix. 


[H]: fluid (pseudo-)stiffness matrix. 
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[M,]: consistent (pSeudo-)fluid mass matrix. 

[L]: coupling matrix at the fluid structure interface (connectivity matrix). 
[O]: zero matrix. 

w: angular frequency. 

p: fluid density. 

c: fluid sound speed. 


T: means transposed. 


ATILA is able to perform: 

‘he Static analysis of an elastic, piezoelectric, or hydroelastic system, 
where the displacement field and/or the electric potential or the pressure fields 
are required. 

Z Modal analysis of an elastic, piezoelectric, hydroelastic system, 
where the eigenfrequencies, the resonance and antiresonance frequencies, and 
the normal modes are computed. 

3. Harmonic analysis of a driven elastic or piezoelectric structure, or 
the scattering of an arbitrary incident wave by an elastic or piezoelectric 
Structure. 

A scattered wave analysis of a spherical pressure wave incident on a solid 
structure is required to investigate the performance of the radiation boundary 
elements, and is presented in the following section. 

2. Harmonic Analysis of a Solid Structure Excited by an 
Impinging Wave 

In this type of analysis, the real and imaginary parts of the pressure field 
(P), the displacement field (U), the potential (®), and the electric current (jaQ) 


are computed. The pressure (P) and the flux pc*¥, are split into incident and 


scattered parts. The normal derivative of the incident pressure on the surface 
boundary S is written 1), where [D] is a matrix obtained by assembling the 


damping elements on the surface, provided by the code. 
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The incident pressure field (P.) is provided via a FORTRAN function 


"INCPRE (x,y,z,k)", included at the end of the main program by the user as 
shown in Appendix B, after [Ref 2]. Utilization of a user-provided incident 
pressure allows the excitation of the structure with an incident wave of any kind, 
while the default function provided with ATILA generates a plane wave e™ 
traveling from the positive to the negative direction along the Ox axis (e™ time 
harmonic dependence is assumed). 

By assigning a specific entry in the ATILA code, we are able to compute 
the total pressure or just the scattered one, utilizing the "TOTAL" or 
"SCATTERED" attribute, respectively. 

The equations of elasticity in the structure, the Helmoltz equation in the 


fluid, and Poisson's equation in the solid material are written in matrix form : 

















aie a eee U F -(LP, 
[Kup]! [K 46] [0] ® |= -Q 
i L [D] oP [H) [MI 
H M res aii , 
pnt ort (Bh Ma) I | erat Gru? on Geta? ~ pte 


where: 


P., : vector of the nodal values of the elastic scattered pressure field, 
P,: vector of the nodal values of the incident pressure field, 


[G]: frequency dependent complex linear matrix, 

[D]: matrix of the damping elements. 

Furthermore, the internal losses of the material used can be taken into 
account throughout a specified program entry "SKYLINE COMPLEX" and when 
no losses are required, with "SKYLINE REAL", respectively. 





Fi 


lll. THREE DIMENSIONAL SPHERICAL MODELS 


A. INTRODUCTION 

In order to evaluate ATILA’s radiation boundary elements, we first 
developed a family of three-dimensional spherical models. These models 
simulate a spherical rigid structure surrounded by an infinite fluid environment. 

A total of five models were employed. All models are composed of an 
inner spherical boundary representing a rigid solid, several concentric spherical 
fluids layers, and an outer spherical boundary composed of radiation boundary 
elements (dampers), each representing a semi-infinite fluid region. 

B. DESCRIPTION OF MODELS 

The models developed are distinguished by their values of the inner 

radius a, the number and thickness L of each fluid layer, and the outer radius R. 


The appropriate dimensionless forms for these parameters are ka, kL, KR, where 
k= is the acoustic wavenumber. Table 1 lists the properties of each model 
Cc 


used in the investigation. 

Model 1 had already been used in the calculation of the transition matrix 
for the scattering of acoustic waves from a thin elastic spherical shell (Ref. 2]. 
This model contains four fluid layers of two different thicknesses. To separately 
investigate the influence of each of the dimensionless parameters (ka, kR, KL), 
four additional models were created. 

Model 2 serves as a reference model. Each of models 3 through 5 is 


distinguished from model 2 in that the value of only one of ka, kL, KR is changed. 






Model! 5 






KR 
1 
0.5 


Number of Fluid Layers 





Mode! 1 Model 2 Model! 3 Model 4 
| cron | ce | ee | 
1 
0 
0 


[Warnunberkinnetes™ | 2 PT 
2 of 
3) 
5 


L 
2 

Radius ain meters 5 1 Ss) 

pamewolAt ISMN iene Chai pasar ee) 
0 1 5 
4 of 3 Of 2 of 2 of 
L=0.5 L=0.5 L= 1.0 L=0.5 
0 


1 
ae, 
enenisseene [0 fs] 


Table 1. Properties of each model used in the investigation. 


e 





0.5 
1.5 





The solid structure is modeled by specifying a zero-flux or boundary 
condition on the scatterer surface (interface between solid and fluid domain). 
This is the default condition. 

For radiation problems, the fluid mesh outer limit must be spherical. 
Therefore, the surrounding solid structure fluid is modeled with a spherical 
surface of dimensionless radius kR = 2.5 or 1.5, for models 2 through 5. 


The fluid (water) is simulated by assuming the following properties: 


il. E=0.222*10'°Pa (Young's modulus) 
D2, p=0.1*10*kg/m* (Density) 
3. v=0.0 (Poisson's Ratio) 


The whole solid structure and fluid mesh was constructed using fifteen- 
node, three-dimensional isoparametric right triangular prismatic elements 
(PRIS15) from the ATILA finite element library [Ref.3]. 

ATILA offers monopolar, dipolar, and multipolar radiation boundary 
elements. Multipolar damping elements were selected to terminate the fluid 


mesh outer radius and the appropriate condition was set in the data input file. 
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The multipolar damping element used is a six-node isoparametric triangular 
element (TRIAO6R) and it has to be spherical and attached to the outer surface 
of the fluid mesh. 

The topology of the prismatic (PRIS15) and of the triangular (TRIAO6R) 
elements are shown in Figures 3 and 4, respectively. The numbers represent 
the nodes and the order of numbering required for each element. 

C. MESH GENERATION 

The major limitation of the three-dimensional finite element model is the 
number of degrees of freedom (D.O.F.) available (displacement along the 
coordinate axes, rotation along the coordinate axes, pressure, electric or 
magnetic potential, and magnetic excitation currents). For this reason the 
models developed were limited by the number of nodes and elements that were 
allowed. 

Specifically, the number of nodes employed in the models ranged from 
1546 to 2346 for the structure and fluid mesh, and the total number of elements 
ranged from 216 to 360. Of these, 144 to 288 are fluid elements and 72 are 
damping elements. The fluid mesh is arranged in successive layers of various 
thicknesses. The original fluid mesh was divided into four layers, two of them 
with 0.25m thickness and the other two fluid layers with 0.5m thickness. 

Figure 5 shows a Mercator projection of the nodes and elements on the 
inner surface of the fluid. The node numbers and selected element numbers 
corresponding to this layer are given. 

For the mesh generation procedure all the nodes were specified using 
spherical coordinates with the origin at the center of the spherical structure. The 
mesh spacing was always less than a quarter of the acoustic wavelength used. 

The aspect ratio of each element according to the reference manual 


should be less than or equal to 3; in our models it is 1 to 4. 


ie 





Figure 3. PRIS15F Element Topology. After Ref.[3] 





Figure 4. TRIAO6R Element Topology. After Ref.[3] 
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Inner Surface Fluid Layer Model, Mercator Projection. 


Figure 5. 
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The internal angles of our elements were modeled between 30 and 90 
degrees while the reference manual states that these should lie between 30 and 
100 degrees. 

The mesh was built in such a way that adjacent elements, like elements 1 
and 2 shown in Figure 5, were able to share common nodes. The top side 
exploded mesh view of the three dimensional model created by ATILA DEPL 
program is shown in Figure 6. 

The types of isoparametric elements used in the mesh generation are 
described in [Ref. 3] and listed in the following Table 2. 


Fluid PRIS15F 15-node triangular 


prism 
Radiation Damping TRIAO6R 6-node triangular 


Table 2. lsoparametric Elements Used in The Three Dimensional Mesh 
Generation 
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Figure 6. Three Dimensional Spherical Model. Topside View 





IV. RESULTS 


A. EXACT ANALYTICAL RESULTS FOR SCATTERED PRESSURE WAVE 
The linearized homogenous wave equation for the propagation of sound 
in ideal (nonviscous) fluids and the time-independent, lossless Helmholtz 


equation for a pressure wave with time-harmonic dependence p(t,r)=p(r)*e, at a 


position r=(r,8,o) and time (t), is given by [Refs. 9, 10]: 


2 
V?p(te)-—5 ED = 0 = VBI) + K7() = 0 (1) 
Cc 


The solution of the Helmoltz Equation (1) for the incident (P,) and 
scattered (P,.) pressure fields in the spherical coordinate system shown in 


Figure 7 is given by Equations 2 and 3: 





Figure 7. The Spherical Coordinates (r,6,0). 


ay 


P(r Oo)= SE Poy her). ¥™(0,4) 
eae eT a n (2) 


CO n 
Pi(.0)= 5 LY Phan hl) (kr)-¥M 6,9) 
n=O0Om=-—n (3) 


where: 


Pint and P,,are the amplitudes of the n,zm components of the incident 


n 
and scattered waves, respectively, 


h (kr) = the nth order spherical Hankel function of the first kind, 
h)(kr) = the nth order spherical Hankel function of the second kind, 
Y™ (6,6) = the spherical harmonic of order n,m, which ts related to the 


associated Legendre polynomial by the equation 


2n+1,(n-m)! 


re Vara -P™M(coso)-el™ , 








Ynm (8, o) zs 


Application of boundary conditions on the inner surface of the spherical 
fluid mesh for vanishing of the pressure or its normal derivative, provides 
solutions for the scattered wave for the case of a pressure release or rigid 


surface. 


1. Analytical Expression for the Scattered Pressure from a 
Pressure Release Spherical Surface 


The following equations apply for the n,m component: 


P._(r,0,0)=¥" (6,9), (kr) (4) 


P..(1,0,0)=[-h (ka) /h') (ka)]- ¥'™(0, 9)! (kr) (9) 
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where: 

a = the radius of the spherical shell, 

kK = the wavenumber of sound, 

i Analytical Expression for the Scattered Pressure from a Rigid 
Spherical Surface 


The following equations apply: 


P._(r,0,9)=¥'(0,9)-h' (kr), (6) 
Pso(t.0.0)=I-(nitkay) /( nf?) (ka) | 1-¥q"(G.0)-MfO (Kn, 1) 


where (h‘(ka))’ and (h!?)(ka))’ are the spatial derivatives of h{ (kr) and h\?)(kr) 
evaluated at the spherical surface. 

Appendix A contains a C program used to calculate analytically the values 
of the pressure scattered/radiated from a vibrating, spherically symmetric 
surface, for the above boundary conditions [Ref. 11]. 


B. NORMALIZED SCATTERED PRESSURE DEVIATION COMPUTED BY 
ATILA 
Using the requried ATILA input data file, the pressure scattered by a rigid 


spherical surface due to an incident spherical harmonic wave was computed. 
The results were compared with an exact analytical solution computed using the 
C program described in Appendix A. Accordingly, the limitations of the boundary 
elements were obtained. 

Since the scattered pressure is proportional to the incident pressure and 
has angular dependence, the error in the ATILA results at each finite-element 
node was quantified by normalizing the deviation from the exact value by the 
maximum magnitude of the scattered pressure at the same radius. This 


represents the normalized scattered pressure deviation computed by ATILA. 
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C. | RESULTS FOR THE THREE DIMENSIONAL SPHERICAL MODEL 
(MODEL 1) 


The following Figures 8 and 9 depict the results of the calculation for the 
original three-dimensional spherical model (model 1) that had been used in the 
previous investigation [Ref. 2]. In this model, the fluid is divided into four layers 
of thickness: 

Ae Layera: 0.25m 

2. Layerb: 0.25m 

oi, Layerc: 0.5 m 

4. Layerd: 0.5m 
The radius of the scattering (inner) surface (a) is 0.5m: the radius of the radiation 


boundary is 2.0m. The frequency (f) is 474Hz; the corresponding wavenumber 


of sound in the water is k=2=2m"! 


Figure 8, from top to bottom, presents the normalized scattered pressure 
deviation computed by ATILA using multipolar dampers for monopolar and 
dipolar incident pressure fields, versus the dimensionless radius (kr) from the 
center of the structure. Figure 9 presents the normalized scattered pressure 
deviation for the quadrupolar incident pressure field. 

The following Table 3 provides the normalized maximum error in 
percentages in the middle of each fluid layer and the angular location of that 
error for the monopolar, dipolar, and quadrupolar incident fields. 

Nodes located at the poles are points with coordinate spherical angles: @ 
(polar) = 000,180 and o (azimuthal) = 000, degrees. Nodes located on the 
equator are points with coordinate spherical angles: 98 (polar) = 090 and o 
(azimuthal) = 000,090,180,270 degrees. Near the equator, nodes are points with 
coordinate spherical angles: 6 (polar) = 075,105 and o (azimuthal) = 
090,180,270 degrees. 
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0.12 
0.1 
0.08 - 


0.06 


__| MONOPOLAR RADIATION 


oe n=0 / m=0 , 
| 


0.02 








0,12 Fiala 
0.1 | 
0.08 } 
0.06 }-— 
0.04 } 
0.02 


DIPOLAR RADIATION 
n=1/m=0 





DIPOLAR RADIATION 
n=1 / m=+1 





Figure 8. Normalized scattered pressure deviation computed by ATILA versus kr 
for monopolar and dipolar incident fields. 
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Figure 9. Normalized scattered pressure deviation computed by ATILA versus kr 
for quadrupolar incident field. 
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Scatterer | Center Center of | Center of |} Center of | Maximum 

Surface of Normalized 

ka=1.0 Layer a =1. =2. =3. Error 
kr=1.25 


Monopolar 2.6 3.8 Or 7.8 5.8 8.6 kr=3.0 | 4.6 
Dipolar bE 2 8.5 11.8 kr=3.0 | 3.9 
n=1, m=0 me | a 


5.8 kr=3.0 | 0.6 


Quadrupolar | 0.3 
ae — E-2: Poles Poles Poles 
Wea: 0.3 0.4 0.5 0.9 Se7 3.6 kr=3.5 0.5 
05 


ee 3.6 kr=3.5 


Dipolar 1.8 PA fo 3.6 &.1 4.0 5.7 kr=3.0 
n=1, m=+1 Equator Equator | Equator Equator Equator Equator Equator 


n=2, m=+2 Near Near Near Near Near Near Near 


Equator Equator | Equator Equator Equator Equator Equator 





Table 3. Percent normalized scattered pressure deviation and location of the 
maximum error for monopolar, dipolar and quadrupolar incident fields. 


From the above table and Figures 8 and 9, it can be concluded that: 


i The greatest normalized errors appear for the following type of 
fields: 
a. - Monopolar, n=m=0: 8.6% at kr=3.0 
b. Dipolar, n=1, m=0: 11.8% at kr=3.0 
C. Quadrupolar, n=2, m=0: 5.8% at kr=3.0 
2. For the above axisymetrical incident pressure fields the 


corresponding location of the maximum error points is close to the poles and 
exactly on the poles for the scatterer’s surface and on the poles for each layer, 


respectively. 
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oo For the cases, where the maximum error points appear on the 
equator or near the equator (dipolar n=1, m=+1, quadrupolar n=2, m=+1), the 
location of the minimum error points is on the poles. 

4. lf the maximum error points are ignored for the monopolar and 
dipolar fields then the maximum normalized error is always less than 4%. 

3 lf the maximum error points are ignored for the quadrupolar field 
then the maximum normalized error is always less than 3.5%. 

6. The normalized scattered pressure deviation computed by ATILA 
as it is presented in Figures 8 and 9 appears to increase when moving further 
away from the acoustic center. 

The above results for the monopolar, dipolar and quadrupolar incident 
fields indicate that, for a given value of n, the maximum normalized errors occur 
for the axially symmetric type of fields (i.e., for m = 0). Hence, in investigating 
the influence of the fluid mesh inner radius (ka), the element's size (kL), and the 
fluid mesh outer radius (KR) on the results, only the axially symmetric incident 
waves were analyzed. 

D. INFLUENCE OF FLUID MESH INNER RADIUS 

In order to evaluate the influence of fluid mesh inner radius on the results, 
models 2 and 3 were used. Recall that model 2 is divided into four equal fluid 
layers of thickness L=0.5m and the scatterer’s radius is 0.5m. Model 3 is divided 
into three equal fluid layers of thickness L=0.5m and scatterer’s radius 1.0m. For 
both cases, the radiation boundary and element's size satisfy kR=2.5 and 
kL=0.5, with k =10m™' and R=2.5m. 

Figure 10 presents the influence of the scatterer’s radius (ka) on the 
normalized scattered pressure deviation, versus kr. The monopolar (a,b), dipolar 
(c,d) and quadrupolar (e,f) fields are shown from top to bottom. On the left side 


the scatterer radius is ka=0.5 and on the right side the scatterer radius is ka=1.0. 
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Figure 10. Normalized scattered pressure deviation computed by ATILA versus 
kr. From top to bottom monopoloar (a,b), dipolar (c,d) and quadrupolar (e,f) 
incident field; on the left: scatterer radius ka = 0.5; on the right: scatterer radius 
ka = 1.0, radiation boundary R and element size L satisfy: kR = 2.5 and kL = 0.5 
for all cases. 
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Table 4 summarizes the results for the maximum error in percent in the 
middle of each layer and the node location for the monopolar, dipolar and 
quadrupolar incident fields when the scatterer radius ka varies and kL and kR do 


not. 







Center of | Center of | Center of | Center of | Maximum 
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Surface 
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Normalized 
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n=m=0 
ka=1.0 
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Dipolar 








3.9 9.8 2.7 
4.9 Ce 0.9 


Poles Poles Poles 


7.0 kr=2.0 
7.8 kr=2.0 


Poles 
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Table 4. Percent normalized scattered pressure deviation and location of the 
maximum error for monopolar, dipolar and quadrupolar incident pressure field 
when ka=1.0 and ka=0.5, radiation boundary kR=2.5, element's size kL=0.5. 


From this table and Figure 10 it is observed that: 
1. For the monopolar field: 
a. The normalized error at a given value of kR decreases as ka 


is increased from 0.5 to 1.0. 


b. The maximum error decreases as the scatterer radius ka is 
increased from 0.5 (7.8%) to 1.0 (7.0%). 

(on For both cases, the maximum error occurs at kKR=2.0. 

d. lf the maximum error points (poles) are disregarded, then 


the maximum error is less than or equal to 4%. 
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2 For the dipolar field: 

a. The normalized error at kr = 1.0 increases by a factor of 5, 
from about 0.3% to about 1.5% as the scatterer radius ka is increased from 0.5 
to 1.0. 

b. The normalized error at all other radii is essentially 
unchanged by varying ka. 

Cc The maximum error still occurs at kr=2.0 as it occurred for 
the monopolar field. 

d. When the maximum error points (poles) are disregarded 
then the maximum error is less than 3%. 

a For the quadrupolar field: 

a. The normalized error at a given node shows very little 
dependence on the value of ka. 

b. The maximum error remains the same (4.1%) as the 
scatterer radius ka is increased from 0.5 to 1.0, from the acoustic center 

Ce For both cases, if the maximum error points are disregarded 
then the error is less than 2%. 

For the above dipolar and quadrupolar fields, the minimum error appears 
on the equatorial points. Moreover, from Figure 10, very similar deviation curves 
for the monopolar (a,b), dipolar (c,d) and quadrupolar (e,f) fields, for ka=0.5 and 
ka=1.0 are observed. Also, we observe that the maximum normalized error on 
the scatterer's surface increases as the scatterer’s radius increases. 

E. INFLUENCE OF FLUID MESH ELEMENT SIZE 

In order to evaluate the influence of fluid mesh element size on the 
results, we used models 2 and 4. Recall that model 2 is divided into four equal 
fluid layers of thickness L=0.5m while model 4 is divided into two equal fluid 
layers of thickness L=1.0m. For both cases, the radiation boundary and 
scatterers radius satisfy kKR=2.5 and ka=0.5, with k=10m™', R=2.5m and 
a=0.5m. 


2/ 


Figure 11 presents the influence of the element’s size (kL) on the 
normalized scattered pressure deviation, versus kr. The monopolar (a,b), dipolar 
(c,d), and quadrupolar (e,f) fields are shown from top to bottom. On the left side 
the element's size is KL=0.5 and on the right side the element's size is KL=1.0. 

Table 5 summarizes the results for the maximum error in percent on the 
scatterer’s surface, in the middle and at the end of each layer and the node 
location for the monopolar, dipolar and quadrupolar fields, when the element’s 


size varies and kR, ka do not. 


Type of scatterer Maximum 
Field Surface Normalized 
ka=0.5 kr=1.0 kr=1.5 kr=2.0 kr=2.5 Error 


Monopolar Poles Poles Poles Poles Poles 


n=m=0 
kL=1.0 4.1 B.S 6.3 4.5 6.3 kr=2.0 
kKL=0.5 3.8 Bal, 7.9 33 7.9 kr=2.0 


Dipolar ~ | Poles Poles Poles Poles Poles 


n=1, m=0 

kL=1.0 1.0 2.6 OES 3.9 3.9 kr=2.5 
kL=0.5 0.5 Pangpe 4.0 3.6 4.0 kr=2.0 
Quadrupolar Poles Poles Poles Poles Poles 
n=2, m=0 

kL=1.0 ; 0.6 0.9 lee 4.2 4.2 kr=2.5 
kL=0.5 ; 0.6 0.7 1.9 4.0 4.0 kr=2.5 





Table 5. Percent normalized scattered pressure deviation and location of the 
maximum error, for monopolar, dipolar and quadrupolar incident pressure field, 
when kL=1.0 and kL=0.5, scatterer radius ka=0.5, radiation boundary kKR=2.5. 
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Figure 11. Normalized scattered pressure deviation computed by ATILA versus 
kr. From top to bottom monopolar (a,b), dipolar (c,d) and quadrupolar (e,f) 
incident field; on the left: element size kL = 0.5: on the right: element size kL = 
1.0. Scatterer radius a and radiation boundary R satisfy: ka = 0.5 and kR = 2.5 
for all cases. 
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From this table and Figure 11 we observe that: 
Ay For the monopolar field: 

a. On the scatterer’s surface ka=0.5, as the element's size is 
increased, the maximum error decreases slightly from 1.7% to 1.2%. 

b. The maximum error at each layer becomes greater as kL is 
increased, except at the maximum deviation point, at kr=2.0. 

‘er For both cases the maximum error occurs at kr=2.0. 

d. Once the maximum error points (poles) are disregarded, 
then the maximum error is less than or equal to 4%. Also, in that case, the 
corresponding deviation for both cases is similar. 

2. For the dipolar field: 

a. On the scatterer’s surface for both cases the deviation is 
almost the same. 

b. As the element's size is increased, the error increases 


except at the maximum deviation point at kr=2.0. 


C. The maximum deviation for kL=0.5 occurs at kr=2.0, while 
for KL=1.0 it occurs at kr=2.5. 
d. When the maximum error points (poles) are disregarded, 


then the maximum deviation is very similar and less than 3% for both cases. 


2. For the quadrupolar field: 


a. On the scatterers surface, ka=0.5, the maximum deviation 
for both cases is 0.2%. . 
b. For both cases, the variation in the normalized deviation 


follows the same trend, and corresponding values are very close to each other. 
on For both cases, when the maximum error points are 
disregarded, the deviation is less than 2%. 
d. The maximum deviation occurs at kr=2.5 for both cases. 
E. INFLUENCE OF FLUID MESH OUTER RADIUS 
In order to evaluate the influence of the mesh outer radius, we used 


models 2 and 5. Model 2 is the reference model and is divided into four equal 
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fluid layers of thickness L=0.5m, while model 5 is divided into two equal fluid 
layers of thickness L=0.5m. For both cases, the scatterer radius a and elements 
size L, satisfy ka=0.5 and kL=0.5, for k = 10m" and a=0.5m. 

Figure 12, presents the influence of the radiation boundary (kR) on the 
normalized scattered pressure deviation versus kr. The monopolar (a,b), dipolar 
(c,d) and quadrupolar (e,f) fields are shown from top to bottom. On the left side 
the radiation boundary is kR=1.5 and on the right side the radiation boundary is 
kKR=2.5. 

Table 6 summarizes the results for the maximum error in percent in the 
middle of each layer and the node location for the monopolar, dipolar and 
quadrupolar fields when the radiation boundary kR varies and kL and ka are the 
same. 

lt is observed from this table and Figure 12 that : 

ue For the monopolar field: 

a. At a given radius, the error increases as the radiation 
boundary kR is increased. For example, in the case where kR=2.5, the error is 
100% greater than the corresponding error when kR=1.5 at the location where 
kr=1.5. 

b. At the outer fluid boundary, the error remains almost the 
same (2%) when the maximum error points (poles) are disregarded. 

C. When kR=1.5, the maximum error occurs at the outer radius 
surface points. On the other hand, when kR=2.5, the maximum error occurs at 
kr=2.0. 

2. For the dipolar field: 

a. The error decreases drastically as we increase the radiation 
boundary kR and becomes almost 100% lower at the boundary for kR=2.5. 

b. The maximum error still occurs at the same radius as it does 


for the monopolar field. 
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Figure 12. Normalized scattered pressure deviation computed by ATILA versus 
kr. From top to bottom monopolar (a,b), dipolar (c,d) and quadrupolar (e,f) 
incident field; on the left: radiation boundary at kR = 1.5; on the right: radiation 


boundary at kR = 2.5; scatterer radius a and element size L satisfy: ka = 0.5 and 
kL = 0.5 for all cases. 
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Scatterer | Center of | Center of | Center of | Center of | Maximum 
Surface Layer d Normalized 
ka=0.5 kr=2.25 Error 
Monopolar Poles Poles 
n=m=0 
kKR=1.5 , 23 2.3 kr=1.5 
kR=2.5 é ; 4.9 ; 7.8 kr=2.0 


Dipolar Poles Poles 


n=1, m=0 

kR=1.5 : 5:3 8.1 kr=1.5 
kR=2.5 . 1.8 , . 4.0 kr=2.0 
Quadrupolar Poles Poles 
n=2, m=0 

kR=1.5 1.6 4.8 kr=1.5 
kKR=2.5 , 0.5 4.1 kr=2.5 





Table 6. Percent normalized scattered pressure deviation and location of the 
maximum error for monopolar, dipolar and quadrupolar incident fields, when 
kKR=1.5 and kR=2.5, scatterer radius ka=0.5 and element size kL=0.5. 

C Once again, if we disregard the maximum error points 
(poles), for KR=2.5, the error is always less or equal to 3%. On the other hand, 
for kKR=1.5, the error is less than 7%. 

d. The minimum error points for both cases appear on the 
equator. 

Sy For the quadrupolar field: 

a. At a given radius, the error decreases as the radiation 
boundary increases. 

b. The maximum error for KR=1.5 and kR=2.5 occurs on the 
outer boundary surface. 

C. For both cases, if the maximum error points (poles) are 


disregarded the error is always less than 2%. 
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d. Finally, the minimum error for both cases appears on the 
equatorial points. 
G. SUMMARY OF RESULTS 

The detailed analyses developed for the five models listed in Chapter III, 
Table 1, indicate that: 

1. For the original three dimensional model (model 1), the maximum 
deviation occurs close to or at the poles for the axially symmetric monopolar 
(8.6%), dipolar (11.8%) and quadrupolar (6%) incident pressure fields. Also, the 
minimum deviation points are at equatorial nodes. Moreover, for the non-axially 
symmetric incident fields, the maximum deviation is less than 5.5% and occurs at 
equatorial points, while the minimum deviation occurs at points located at the 
poles. 

2. Further investigation of the influence of the fluid mesh inner radius 
(ka), the element's size (kL), and the fluid mesh outer radius (KR), on the results 
for axially symmetric incident pressure fields, revealed that: 

a. For monopolar and dipolar incident fields, when the poles 
are disregarded: 

(1) As the radiation boundary kR increases while ka and 
kL remain constant, the maximum deviation decreases by approximately the 
Same amount in percent. This occurs specifically from 6.6% when kR=1.5, to 
4% when kR=2.5, for ka=0.5 and kL=0.5, respectively, which corresponds to 
65%. 

(2) As the scatterer’s radius ka increases, while kL and 
kR remain constant, the maximum deviation decreases by a small amount in 
percent. This occurs, specifically, from 4% when ka=0.5, to 3.7% when ka=1.0, 
for KL=0.5 and kR=2.5, respectively, which corresponds to 8%. 

(3) As the element’s size kL increases, while ka and kR 
remain constant, the maximum deviation remains constant. 

b. For a quadrupolar incident field, the maximum deviation 


remains essentially the same between 1.9% and 1.8%. 
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Table 7, which combines three tables, summarizes our results for the 
maximum normalized error in percent, when the poles are disregarded, and the 
type of field for which this error occurs. Also, from this table it can be observed 
that the maximum deviation occurs for monopolar and dipolar axially symmetric 


incident fields. 





ka=0.5 
4.0 4.0 4.0 a7 ka=1.0 Sle 
ET 


Table 7. Percent normalized maximum scattered pressure deviation and type of 
incident pressure field, when poles are disregarded, computed by ATILA. 


ka=0.5 kL=0.5 kL=1.0 kL=0.5 ka=0.5 ka=1.0 
i i | ee 
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V. CONCLUSIONS AND SUGGESTIONS FOR FURTHER 
INVESTIGATION 


A. CONCLUSIONS 

Five three-dimensional spherical models were developed in order to 
evaluate the radiation boundary elements used in the ATILA finite element code. 
The models simulate a rigid spherical solid structure surrounded by an infinite 
fluid. 

Monopolar, dipolar, and quadrupolar incident spherical pressure wave 
fields were imposed, and the scattered pressure was calculated using ATILA. 
Since the scattered pressure is proportional to the incident pressure and has 
angular dependence, the error in the ATILA results at each finite-element node 
was quantified by normalizing the deviation from the exact value by the 
maximum magnitude of the scattered pressure at the same radius. This 
represents the normalized scattered pressure deviation computed by ATILA. 

The range of values of the dimensionless lengths which characterized the 
problem was: 

The fluid mesh inner radius ka (scatterer radius): 2.0 , 1.0 , 0.5 
The fluid mesh element size kL: 0.5 , 1.0 
The radiation boundary KR: 4.0 , 2.5, 1.5 

The maximum deviations observed were for the axially symmetric 
monopolar (9%), dipolar (12%), and quadrupolar (6%) incident pressure fields. 
For these cases, it can be concluded that the maximum deviation points are at 
the poles, while the minimum deviation occurs on the equatorial points. 

When the poles are disregarded, for the monopolar and dipolar incident 
pressure fields there is a strong influence of the fluid mesh outer radius (KR) and 
a weak influence of the fluid mesh inner radius (ka) on the results. Specifically, 
as the fluid mesh outer radius increases and/or the fluid mesh inner radius 


increases then the normalized scattered pressure deviation decreases. For the 


3/ 


quadrupolar field, the maximum deviation remains essentially constant, 
independent of the above factors (ka, kR). Furthermore, variation of the 
elements size (kL) does not affect the maximum normalized scattered pressure 
deviation. 

Moreover, for the non-axially symmetric incident pressure fields, the 
maximum deviation observed was less than 5.5%, and was found to occur at 
nodes on the equator, while the minimum deviation points were located at the 
poles. 

B. SUGGESTIONS FOR FURTHER INVESTIGATION 

The magnitudes of the errors in the scattered pressures computed by 
ATILA were larger than expected, based upon the mesh sizes employed in the 
models; for all of the models the radial length L of the fluid elements was less 
than 1/6 wavelength. It would have been desirable to examine the calculated 
scattering for a more refined mesh. However, to divide the element scale size by 
a factor of 2 would have resulted in models which exceed the allowed number of 
degrees of freedom. Of course, this is because all the models employed were 
fully three-dimensional, and did not make use of any possible reductions in size 
due to wave function symmetry (in fact, our choice of e” for the azimuthal wave 
function precluded reduction in that direction). Also, we are interested in the 
computed scattering in all of the spherical harmonic wave components for an 
incident wave of only one component. 

It is suggested that further investigation of the performance of the ATILA 
radiation boundary elements be conducted using two-dimensional models. _ It 
might also be advantageous to use sin @ or cos 4 instead of e'"’ to represent the 
azimuthal dependence. This will allow a finer mesh to be employed, and a 


broader range of values of ka, kL, and kR to be examined. 
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APPENDIX A. C CODE FOR ANALYTICAL PRESSURE 
SCATTERED CALCULATION 


KKEKKEKEKEKKEEKEKEKKEKK KEKE KEKE KKK KEKE KEKE KEKRER KERRI 


* PROGRAM SCHAR.C , by Panagiotis Sinanoglou 2/2/1996 
* Program:Computes the Scattered Spherical Harmonic pressure from a rigid 
* spherical structure at the point r1(x,y,z), for the Wavenumber k, for: 
* natural boundary and pressure release condition 
* Input Variables: 
* x,y,z: Cartesian Coordinates of the Node Point. 
k: Wavenumber. 
r2: Radius of the Scatterer's Surface. 
n,m: Orders of Spherical Hankel and Legendre functions. 
Output Variables: 
r1: Radius in meters. 
phi: Azimouthial angle in degrees. 
theta: Polar angle in degrees. 
rp5: Real part of computed scattered pressure at r1 in pascals 
rp6: Imaginary part of computed scattered pressure atr1 in pascals 
* Formula for Incident Wave: Pinc=Pnm(costheta)*e’(imphi)*hn(1)(kr2) 
*Formula for Scattered Wave from Rigid Boundary: 
* Psc=-Pnm(costheta)*e“(imphi)*[hn'(1)(kr2)/hn'(2)(kr2)}*hn2(kr1) 
* Formula for Scattered Wave from Pressure Release Surface: 
* Psc=-Pnm(costheta)*e“(imphi)*[hn(1)(kr2)/hn(2)(kr2)]*hn2(kr1 ) 
* External Functions: sphbes.c, plgndr.c, bessjy.c, beschb.c, chebev.c, 
complex.c, nrutil.c 
* Header Files for prototype function declaration:nr.h nrutil.h,complex1.h 
* Functions: 
* Normalized Legendre[ plgndr(n,m,x) ],Spherical Hankel for the 
* real and the imaginary parts[ sphbes(n, (k*r2),&xsj,&xsy,&xsjp,&xsyp) J, 
* Complex(structures)[ Cexp(mphi), Cdiv(a,b), Cmul(a,b), RCmul(a,b) ] 
* File nphra33n.c provides node coordinates(x,y,Z). 
* File dim1.dt stores the corresponding spherical coordinates(r1,theta, phi) 
* File dim13.dt stores the Legendre plgndr(n,m,(z/r1). 
* File dim11.dt stores the real,imaginary part of exp. function e“(imphi). 
* File dim14.dt stores the ratio of the derivatives of the Spherical Hankel for the 
real and 
* imaginary part of hn'(1)(kr2)/hn'(2)(kr2) for natural boundary conditions and the 
ratio of 
*hn(1)(kr2)/hn(2)(kr2) for pressure release,evaluated at the scattering surface, at 
r2=0.5m. 
* File dim15.dt stores the real and the imaginary part of the Spherical 
* Hankel hn2(kr1) at the node range r1. 


os) 


+ + + + + + + + 


* 


* File dim16.dt stores the real and the imaginary part of the scattered pressure 
Psc, at the * range r1. 


KRKEEKKEKKKKKEKKKEKKEEKEKKKEKKEKKEKEKKREEKKKEKKKKKKKKK KKK KKK KKKKKKKKKKKKKKEKK 


#include <stdio.h> 
#include <math.h> 
#include<string.h> 
#include<stdlib.h> 
#include “complex1.h" 
#define NRANSI 
#include “nr.h" 
#include “nrutil.h" 
#define pi 3.141592654 


main() 


Alle “hs Gili iS: .Gee. Zee PIE oi4 BILE *f5: Pile siGslee sa. 
FILE *f20; FILE *f21;FILE *f22: 

fcomplex a,b,hn2; 

float r2=0.5,k=2.0; 

float sj,Sy,Sjp,Syp,xsj,xsy,xsjp,xsyp,rp1,rp2,rp3,rp4,rp5,rp6; 

int j,rpm; 

float fac, val,t,f11,f12,f13,f14,Pnm,s1,s2: 
unsigned int factorial(unsigned int a); 
double x,y,z,Z1,r,r1,theta,phi,g,mphi:; 

int d,i,N,ch,m,n; 

char g[11],h[11],c[11],e[1 1], f[1 1]; 


f3=fopen("dim1.dt","w"); 
printf("Enter the integer valuesfor N,n,m\n"); 
scanf("%d %d %d",&N,&n,&m); 
f4=fopen("sphra33b.c","r"); 
f5=fopen("dim11.dt","w’); 
f6=fopen("dim12.dt","w’'); 
f7=fopen("dim13.dt","w"); 
f20=fopen("dim14.dt","w'); 
f21=fopen("dim15.dt","w"): 
f22=fopen("dim16.dt","w"); 
for(i=O; I<N; ++i) 


fscanf(f4 ,"Yolf Ylf “lf Yd 
~AlfYshshHsVs*Ms in" &xX, &y,&Z,&d,&g,&g,&h, &c, &e, &f); 
/* Converts from cartesian(x,y,z) to spherical coordinates(r1 ,theta,phi).*/ 
if((x>0.0 && y>=0.0)) 
{ r1=saqrt(x*xty*y+z*z); 
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phi=atan(y/x); 
theta=acos(z/r1):} 
else if(x==0.0 && y==0.0 && z==0.0) 
{r1=0.0; theta=0.0; phi=0.0;} 
else if (x==0.0 && y==0.0 && z!=0.0) 
{ r1=sqrt(x*xt+y*y+z*z); 
phi=0.0; 
r1=sqrt(x*xt+y*y+z*z); 
theta=acos(z/r1);} 
else if(x==0.0 && y>0.0 &&(Z==0.0]|z!=0.0)) 
{ phi=pi/2; 
r1=sqrt(x*xt+y*y+z*z); 
theta=acos(z/r1);} 
else if(y==0.0 && x<0.0 ) 
{ phi=pi; 
r1=sart(x*xty*y+z*z); 
theta=acos(z/r1);} 
else if(y<0.0 && x==0.0) 
{ phi=3*(pi/2); 
r1=sart(x*xty*y+z*z); 
theta=acos(z/r1);} 
else if ((x<0.0 && y>0.0) ) 
{phi=+(pi)+atan(y/x); 
r1=saqrt(x*xt+y*y+z*z); 
theta=acos(z/r1); } 
else if (x<0.0 && y<0.0) 
{ phi=pitatan(y/x); 
r1=sqrt(x*xt+y*y+z*z); 
theta=acos(z/r1);} 
else if (x>0.0 && y<0.0) 
{ phi=(2*pi)t+atan(y/x); 
r1=sqrt(x*xt+y*y+z*z); 
theta=acos(z/r1);} 


theta=theta*(180/pi); 
phi=phi*360/(2*pi); 
fprintf(f3," %+4.3If %+6.3If %+6.3lf\n",r1,(theta),(phi)); 


/* The required value of xx=cos(theta) for the Pnm(xx) Is: (z/r1) */ 
/* Calculates the Normalized Legendre function “pigndr(n,m,(z/r1) */ 
fprintf(f7,"%4s %4s %10s %24s\n","n","m","x","plgndr(n,m,x)"); 
rom=abs(m); 
if(m>=0 ) 
{ f13=factorial(n-m), 
f14=factorial(n+m); 
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$1=sqrt((2*n+1)*f13/(4*pi*f14)); 
Pnm=(plgndr(n,m,(z/r1))*s1): 
fprintf(f7 ,"%4d %4d %13.6If %19.6e\n",n,m,(z/r1),Pnm); 


} 

if(m<0Q) 

{ f13=factorial(n-rpm); 
f14=factorial(nt+rpm); 
s2=sqrt((2*n+1)*f14/(4*pi*f13)): 
f11=factorial(n-rpm); 
f12=factorial(n+rpm); 
r=pow(-1,rpm); 

t=r*(f11/f12); 
Pnm=(plgndr(n,rpm,(z/r1))*t*s2); 

fprintf(f7 ,"%4d %4d %13.6lf %19.6e\n",n,m,(z/r1),Pnm); 


/* Calculates the e“(imphi) Real/Imagin parts */ 


mphi=(m*phi)*2*(pi/360); 

fprintf(f6,"%lf lf \n",Cexp(mphi).r,Cexp(mphi).i); 
/* Calculates the (z/r1),phi(degrees) */ 

fprintf(f5,"%+6.3lf %+6.3lf \n",(z/r1),(phi)): 
/* Calculates the ratio of the Spherical Hankel for the real and the imaginary 
parts */ 
/* of hn'(1)(kr2)/hn'(2)(kr2) on the scattering surface at r2=0.5 for N.B.C , or the 
a 


/* ratio of hn(1)(kr2)/hn(2)(kr2) for pressure release condition 
*/ 


sphbes(n,(k*r2),&xsj,&xsy,&xsjp, &xsyp); 


a.r=xsjp; /* a.r=xsj; */ 
a.l=xsyp; /* a.l=xsy; */ 
b.r=xsjp; /* b.r=xsj; */ 
b.i=-xsyp; /* b.i=-xsy; */ 


fprintf(f20 ,"%30s \n","The real and the imag. of hn'(1)/hn‘(2) is:"); 
fprintf(f20,"“%f %f \n",Cdiv(a,b).r,Cdiv(a,b).1); 
/* Calculates The real and the imag. of hn2(kr1) */ 
sphbes(n,(k*r1),&xsj,&xsy,&xsjp,&xsyp); 
hn2.r=xsj; 
hn2.i=-xsy; 
fprintf(f21,"%30s \n","The real and the imag. of hn2(r1) is:"); 
fprintf(f21,"%f %f\n",hn2.r,hn2.1); 
/* Calculates The real and the imag.parts of the scattered pressure at the */ 
/* radius (r1) :Psc=-Pnm(costheta)*e“(imphi)*[hn'(1)(kr2)/hn'(2)(kr2)]*hn2(kr1)*/ 
rpo5=Cmul(RCmul(-Pnm,Cmul(Cdiv(a,b),Cexp(mphi))),hn2).r; 
rop6=Cmul(RCmul(-Pnm,Cmul(Cdiv(a,b),Cexp(mphi))),An2).1; 
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fprintf(f22,"%f Yf\n",rp5,rp6): 


fclose(f4); 
fclose(f3); 
fclose(f5); 
fclose(f6); 
fclose(f7); 
fclose(f20); 
fclose(f21); 
fclose(f22); 
return 0; 


} 
#undef NRANSI 


unsigned int factorial(unsigned int a) 


if((a==1)||(a==0)) 
return 1: 
else 
{a*=factorial(a-1); 
return a; 


} 


KREKKKKKKEKREKKEKEKKKKKEEEEKKKKEEEKEKEKKKEEKKKEEEKEKRKKKKKKKKKKK KKK KKK 
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APPENDIX B. FUNCTION INCPRE(X,Y,Z,K) 


FUNCTION INCPRE(X,Y,Z,K) 

* PROGRAM BY ARTHUR LOBO DA COSTA RUIZ 12/23/93 

* MODIFIED BY STEVEN R. BAKER 3/4/96 

* FUNCTION: 

* COMPUTES THE INCIDENT SPHERICAL HARMONIC PRESSURE AT THE 
“POINT (X,Y,Z),FOR THE WAVENUMBER K 

* VARIABLES INPUT: 

X,Y,Z: CARTESIAN COORDINATES OF THE POINT 
d K: WAVENUMBER 

* VARIABLES OUTPUT: 

RADIUS: R INMETERS 

PHI INDEGREES (HORIZONTAL ANGLE) 

THETA IN RADIANS (AZIMUTAL ANGLE) 

INCPRE : PRESSURE AT POINT IN PASCAL 


FORMULA USED: 


INCPPRE= H1(N)*P(M,N)*DCMPLX(DCOSD(M‘PHI),DSIND(M*PHI)) 
THIS MAKES AN OUTGOING INCIDENT WAVE FOR E*IWT TIME 
DEPENDENCE , AS FOR ATILA 
DOUBLE PRECISION K,X,Y,Z,R,PHI, THETA,KR 
REAL*8 F1,PMN(-2:2,0:2),_LOUT 
INTEGER N,M,NMAX 
COMPLEX*16 INCPRE,H1(0:2),H1OUT 
*N AND M ARE ORDERS FOR HANKEL AND LEGENDRE FUNCTIONS 
N=0 
M=0 
* H1OUT AND LOUT ARE HANKEL AND LEGENDRE OUTPUTS 
* REF TONANDM 
NMAX=2 
* TRANSFORM CARTESIAN COORDINATES (X.Y,Z) 
* INTO SPHERICAL COORDINATES 
* R(RADIUS), PHI(AZIMUTAL ANGLE) AND THETA(POLAR ANGLE) 
R=DSQRT(X*X+Y*Y+Z*Z) 
PHI=DATAN2D(Y,X) 
IF((X.EQ.0).AND.(Y.EQ.0)) PHI=0.0D0 
THETA=DACOS(Z/R) 


+ + + + + + + HH 4 + 


45 


KR=K*R 

* NMAX IS THE MAXIMUM NUMBER OF HARMONICS 
CALL HANKEL1(KR,NMAX,H1) 

* SUBROUTINE HANKEL1 RETURNS SPHERICAL HANKEL 

* FUNCTIONS OF THE FIRST KIND JN +1 YN 
CALL LEGNDR(THETA,NMAX,PMN) 
LOUT=F1(N,M)*PMN(M,N) 

* SUBROUTINE LEGNDR RETURNS ASSOCIATE LEGENDRE FUNTION 
H1OUT=H1(N) 
INCPRE=H10UT*LOUT*DCMPLX(DCOSD(M*PHI), DSIND(M*PHI)) 
IF ((R.LE.0.501).AND.(R.GT.0.499)) THEN 

PRINT *,"X,Y,Z =",X,Y,Z 

PRINT *,PHI, THETA 

PRINT *,"K,R,KR =",K,R,KR 

PRINT *,"REAL HANKEL1 =",H10UT 

PRINT *,"LEGENDRE =",LOUT 

PRINT *,INCPRE 

PRINT se ia a halal 
ELSE 
CONTINUE 
ENDIF 
RETURN 
END 

C KREKEKKKKKKKKKKEEKKKEKKKKKKEKKEEKEEKEKKEKEEEEREKKKKKKEEKEKEKEEREEKEERKKKEKKEEEK 

SUBROUTINE HANKEL1(X,NMAX,H1) 
IMPLICIT REAL*8 (A-H,O-Z) 
COMPLEX*16 H1(0:NMAX) 


+ * + + %+£ ££ & 


C GIVEN THE VARIABLE X, AND THE MAXIMUM ORDER NMAX, 
C THIS ROUTINE GENERATES THE SPHERICAL HANKEL FUNCTION OF 
THE 
C FIRST KIND HiN FOR ALL N FROM 0 TO NMAX (INCLUSIVE) 
C INPUT: 
C X=DOUBLE PREC. VARIABLE (RADIUS) 
C NMAX = INTEGER MAXIMUM ORDER OF BESSEL FUNCTIONS 
DESIRED 
OUneUN: 
H1(N) = ARRAY OF SPHERICAL HANKEL FUNCTIONS H1N(X), WHERE 
HiN=JN+I1YN 
THIS ROUTINE IS BASED ON THE RECURSION FORMULA 
FROM ABRAMOWITZ & STEGUN: 10.1.10 & 10.1.15, PP.438-9 
THE F'S ARE THE COEFFICIENTS OF ORDER N & -(N+1), 
THE FO'S ARE OLD F'S, FOR RECURSION 
IF (X .LE. 0.0D0 ) THEN 


70) OO Ge @ 
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H1(0) = DCMPLX(1.0D0,-1.0D35) 
DO 2 N= 1, NMAX 
H1(N) = CMPLX(0.0D0,-1.0D35) 
CONTINUE 
RETURN 
END IF 
SX = DSIN(X) 
CX = DCOS(X) 
XINV = 1.0D0/X 
M1N = -1.0D0 
FN = XINV 
FMN = 0.0D0 
FNO = FMN 
FMNO = FN 
DO 4. N=0, NMAX 
H1(N) = CMPLX( FN*SX + M1N*FMN*CX, -FN*CX + M1N*FMN*SxX ) 
T1 = (2*N+1)*XINV 
T2 = T1*FN - FNO 
FNO =FN 
FN =T2 
T2 = -T1*FMN - FMNO 
FMNO = FMN 
FMN = T2 
M1N =-M1N 
CONTINUE 
RETURN 
END 


C KRKEKEKKEKKEEREEKREKEEEREREREREKEKEREREEREEREEEKREKEEEEEEEEEEEREREEE 


SUBROUTINE LEGNDR(THETA,NMAX,PMN) 
IMPLICIT REAL*8 (A-H,O-Z) 
REAL*8 F1,PMN(-NMAX:NMAX,0:NMAX) 


GIVEN THE VARIABLE THETA, AND THE MAXIMUM ORDER NMAX, 
THIS ROUTINE GENERATES THE ASSOC. LEGENDRE FUNCTIONS PMN 
OF THE ARGUMENT COS(THETA) (THETA MUST BE BETWEEN 0 & PI) 
FOR ALL N FROM 0 TO NMAX (INCLUSIVE) 
AND FOR ALL M FROM -N TO N (SOME OTHERS SET TO ZERO) 
INPUT: 

THETA = VARIABLE (POLAR ANGLE), MUST BE BETWEEN 0 & PI 


C NMAX= INTEGER MAXIMUM ORDER OF LEGENDRE FUNCTIONS 

C DESIRED 

C OUTPUT: 

C PMN = DOUBLE PREC. ARRAY, CONTAINS ASSOC. LEGENDRE FNS 
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C THIS ROUTINE IS BASED ON THE RECURSION FORMULAE 
C FROMABRAMOWITZ & STEGUN 
X = DCOS(THETA) 
SINTHT = DSIN(THETA) 
IF (SINTHT .GT. 0. ) THEN 
SININV = 1.0D0/SINTHT 
ELSE 
SININV = 0.0D0 
END IF 
C SET VALUES FOR N =0, 1 (NMAX MUST BE AT LEAST 1) 
PMN(0,0) = 1.0D0 
PMN(1,0) = 0.0D0 
PMN(-1,0) = 0.0D0 
PMN(0,1) = X 
PMN(1,1) = -SINTHT 
PMN(-1,1) = SINTHT*0.5D0 
C IN LOOP, TNP1 = 2*N+1, TNP2FC = (2*N+2)!, M1N = (-1)**(N+1) 
TNP1 = 1.0D0 
TNP2FC = 2.0D0 
M1N = -1.0D0 
DO 4.N = 1, NMAX-1 
TNP1 = TNP1 + 2.0D0 
TNP2FC = TNP2FC * TNP1 * (TNP1+1) 
M1N =-M1N 
DO 3M=-N,N 
PMN(M,N+1) = (TNP1*X*PMN(M,N) - (N+M)*PMN(M,N-1))/(N-M+1) 
3. CONTINUE 
PMN(N+1,N) = 0.0D0 
PMN(-N-1,N) = 0.0D0 
PMN(N+1,N+1) = (X*PMN(N,N+1) - TNP1*PMN(N,N)) * SININV 
PMN(-N-1,N+1) = M1N*PMN(N+1,N+1)/TNP2FC 
4 CONTINUE 
RETURN 
END 
FUNCTION F1(NN,MM) 
IMPLICIT REAL*8 (A-H,O-Z) 
REAL*8 FACT1,FACT2,PI 
P1=4,0D0*DATAN(1.0D0) 
FACT1=1.0D0 
DO 10 1=1,NN+MM 
10 FACT1=FACT1°*! 
FACT2=1.0D0 
DO 20 1=1,NN-MM 
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20 FACT2=FACT2*I 
F1= DSQRT((2*NN+1)*FACT2/4.0D0/PI/FACT1) 
RETURN 
END 
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